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Abstract. 

The stability of a low temperature Bose-Einstein condensate with attractive 
interactions in one and three dimensional double-well potentials is discussed. In 
particular, the tunnelling dynamics of a condensate under the influence of a time- 
dependent potential gradient is investigated. The condensate is shown to collapse at 
a critical potential gradient which corresponds to a critical number of atoms in one of 
the two wells. The sensitivity of this tunnelling induced collapse could provide a useful 
tool in the study of condensates with attractive interactions. 
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1. Introduction 

The experimental realization of Bose-Einstein condensation (BEC) in dilute atomic 
gases of 87 Rb [I], 23 Na 0, 7 Li 0, X H gj, metastable 4 He 0, 85 Rb 0, 41 K [7j, 133 Cs 
jHj and 174 Yb jHj has stimulated theoretical studies of the properties of ultracold Bose 
gases. The densities of such systems are typically sufficiently low, that their interactions 
can be described by a single parameter, the s-wave scattering length a. In addition 
one can control not only the strength of these interactions but also whether they are 
attractive (a < e.g. 7 Li, 85 Rb, 133 Cs at low magnetic fields) or repulsive (a > 0) 
using magnetic-field induced Feshbach resonances ^Oj • Whereas a harmonically confined 
BEC with repulsive interactions is stable for any number of atoms, a condensate with 
attractive interactions is only stable if the atom number M is smaller than a critical 
value M a EED EE1 ECU El US EEJ3 EE2 CHI. For M > M a the interaction energy 
exceeds the zero-point kinetic energy and a collapse occurs, as extensively studied in 
[IHI EDI EH E2I ESI EH E3 EEl E3 EH| In experiments on 7 Li, the magnetic field is held 
fixed and the number of condensate atoms grows up to M cr where a partial collapse 
occurs jinilHni- In this system, it should be noted that the attractive interactions can 
lead to the formation of bright matter-wave soliton trains in elongated optical traps 
[3*T] . In contrast to the 7 Li experiments, in the case of 85 Rb, a Feshbach resonance 
[Oj has been used to switch the scattering length from positive to negative values, 
producing a condensate with M M cr which subsequently collapses [22], as modelled 

in [SHI ES ESI EF3 EH EHl EH] 

Although much work has been done to study the collapse properties of dilute BECs 
in a single-harmonic trap, an interesting subject is the behaviour of a condensate with 
attractive interactions in a double-well potential. Magnetic [HJ EE] and optical [JHJ 
double-well potentials have been created in recent experiments and a proposal for a 
magnetic double- well has been reported [45J. Although condensates with a > in a 
double-well system have received considerable theoretical attention jlHl 113 HH1 HH EUJ 
IHU El] ) , only Adhikari jHH] and Coullet et al. [HE] have studied the a < case. 

In this paper, we investigate the stability of a low temperature atomic condensate 
with attractive interactions in a double- well potential, where the Josephson dynamics are 
induced by the addition of a time- dependent potential gradient. Our analysis is based 
on the Gross-Pitaevskii equation. The aim is to discuss the stability of the condensate 
in the double-well potential and show how tunnelling between the wells can lead to 
collapse. We do not, however, intend to describe the collapse dynamics, for which, one 
should at least additionally incorporate a suitable 3-body loss term (3S1 HO] ■ 

Section 2 introduces the main formalism and briefly reviews the eigenenergies of 
the double-well potential as a function of the potential gradient in both one and three 
dimensions. Section 3 discusses the stability of a BEC in a three dimensional symmetric 
double- well as a function of the barrier height. Finally, Section 4 considers the time- 
dependent evolution of a 3D system as the gradient is increased at a constant rate. We 
find that a collapse occurs as the gradient is increased above a critical value, which, 
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within our formalism, can be predicted by the eigenstate curves. 



2. Eigenenergy levels of a BEC in a double-well potential 

At low temperatures, the condensate wavefunction ip (r,t) (normalized to unity) obeys 
the dimensionless Gross-Pitaevskii (GP) equation, 



■ d , , , 



±V 2 + V(r)+g m \ij(r,t)f 



Hr,t). (1) 



In the above equation, V (r) is the confining potential and c^d = g/(a? L nuj±), is the 
dimensionless effective interaction term, where g = J\f(4:iih 2 a/m) is the usual three- 
dimensional scattering amplitude, defined in terms of the s-wave scattering length a, 
the total number M of atoms of mass m, and the harmonic oscillator length in the 



transverse direction(s) a± = yJh/moj±, with oo± the corresponding trapping frequency. 
The confining potential, 

V (r) = ~ (p 2 + A V) + h exp (-z 2 ) + 8z . (2) 

describes an axially symmetric potential, with p 2 = x 2 + y 2 and asymmetry parameter 
A = uj z /uj± plus a Gaussian potential of height h, located at z = 0. In addition, a linear 
potential 5z of gradient 5 is pivoted at the centre of the trap. For 5 > 0, the right well 
has higher potential energy and the trap centre is additionally shifted into the z > 
region; however, this shift is negligible for the parameters studied throughout this work. 
Throughout this paper we work in dimensionless (harmonic oscillator) units, by applying 
the following scalings: space coordinates transform according to r\ = a±~ 1 r i , time 
t' = ui±t, condensate wavefunction ip' (r', t') = \fa^ip (r, t) and energy E' = (huj ± Y L E 
(primes henceforth neglected for convenience). 

We can also consider the ID form of Eq. (JT|), which accurately takes into account 
the transverse dynamics of 'cigar' condensates (A < 1), 



i>{z,t), (3) 



where ip(z) is normalized to unity, gin = ^30/(2^0^) is the ID self-interaction parameter 
that matches the ID and 3D axial density profiles and Vid(z) is the confining potential 
in the axial direction. 

V m (z) = \z 2 + h exp (-z 2 ) + 5z . (4) 

The eigenstates of the double-well condensate in ID and 3D are calculated by 
substituting ip (r, t) = e~ lfJ ' t ^ (r), where \x is the chemical potential (of the one and three 
dimensional system respectively), and solving the resulting time-independent equation 
as discussed in j^HI- For small nonlinearities there are only two levels and, for 5 = 0, 
the eigenstates are a symmetric ground state ty g and an antisymmetric first excited 
state, \l/ e , with equal population in both wells. Sufficiently large interactions lead to 
the appearance of a loop structure (see e.g. [301 157]). This loop structure appears in 
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Figure 1. Eigenenergies /i (fko±) for the double- well for an attractive Bose gas as a 
function of the potential gradient <5 (huj±/a±) indicating the self-interaction energy, Ec, 
and the Josephson coupling energy, Ej in each case. The horizontal dotted grey line 
corresponds to the zero energy of the two-state model. The value of the nonlinearity 
used here is #3D = — tt, corresponding to Eq — — 0.379hu± and Ej = 0.071Thlj±. The 
eigenstates at the centre of the trap are also shown in each case. We assume a spherical 
trap geometry (A = 1), and a Gaussian barrier of height h = 4hui± located at the centre 
of the trap. 



the first excited state for a > 0, and in the ground state for a < 0. In this paper we 
consider the latter case (Fig. 1), for which the corresponding wavefunctions at 5 = 
are: (i) a symmetric ground state with equal population in both wells, (ii) an anti- 
symmetric state with equal population in both wells and a phase difference of tc across 
the trap centre, which we shall henceforth refer to as \l/ e , and (iii) two lower energy, non- 
symmetric ground state solutions with most of the particles in either the left (dashed) 
or the right (dotted) well (the so-called 'self-trapped' states jUj EBl E3 El ) • 

Under certain conditions, the eigenenergies \i can also be reproduced by the two- 
state model (see e.g. [IBJ HH EB1 EHl E0j) described by the Hamiltonian, 

H= \{ '^-E^ A - BcJV ) ' 

where N = [Ng — N r )/Af is the fractional population difference, A is the potential 
energy difference between the left (£) and right (r) wells, Eq = g(^e, r \ \^e,r\ 2 \^e,r) is the 
self-interaction energy, g is the nonlinearity and Ej = — 2(^| (— |V 2 + ^5=0) \^ r ) is the 
Josephson coupling energy. The energy splittings E c and Ej are indicated in Fig. 1. For 
the case of a potential gradient 5, we can write A = a5, where a is a numerical factor 
determined numerically from the GP solution. The picture for the eigenenergy levels as 
a function of 5 shown in Fig. 1, is only valid in ID and 3D for a range of nonlinearities. 

Focusing first on the ID case, we plot in Fig. 2(a) the two-state model parameters 
I .Eel and Ej at 5 = (symmetric double- well) as a function of gm- For #id < 0, 
increasing the magnitude of the nonlinearity leads to the appearance of a loop structure 
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Figure 2. (a) The self-interaction energy \Ec\ (clashed lines) and the Joscphson 
energy Ej (solid) at S — as a function of the nonlinearity gi^ for the ID confining 
trap with h = Ahuj± . (b) Typical eigenenergies /j, as a function of the potential gradient 
6 for various nonlinearities within regions I-III. Shown are (from left to right) the cases 
giD = -4, -0.5, and 0.5. For g 1D < -3.79, Ej < 0. The horizontal dotted grey line 
corresponds to E = 0. 



in the ground state at a critical point \Eq\ = Ej. As a ID condensate is stable against 
collapse [3J [TTJ Uni I23 UH I2H Ell, l-^cl remains finite and a loop structure is always 
observed. Note that, for a nonlinearity less than a critical value, the splitting Ej 
becomes negative, signifying an inversion of the lowest two energy eigenvalues. This 
can be explained using // = E - \/2\g\ J \^\ 4 dr = E - |£ int |. Although E is larger 
for the first excited state than for the ground state, the first excited state has more 
negative interaction energy (as its peak density is higher), thereby reducing [i below 
the ground state value. Fig. 2(b) shows typical eigenenergy levels for the ground and 
the first excited states as a function of the potential gradient for different values of the 
nonlinearity. The two-state model correctly reproduces the eigenenergy curves for small 
values of \g\v>\ but cannot reproduce the inversion of the eigenenergy levels for g 1D <^ 1. 
The case gin > has been discussed in our earlier work [58] . 

In Fig. 3 we consider the 3D case. In Fig. 3(a) we plot the splittings \Eq\ and Ej 
at 5 = as a function of g-^. In contrast to ID, for g 3 rj < the condensate collapses 
when the atom number, or magnitude of the nonlinearity, exceeds a critical value. The 
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Figure 3. (a) As in Fig. 2, but for the 3D case, with A = 1 and h = 4hu>±. (b) 
Typical eigenenergies shown here are for each of the four regions FIV in (a), namely 
(from left to right) (?3D = — 6, — 7r, and it. Compared to the ID case (Fig. 2), in 
3D there is an additional region (I) corresponding to the case where the self-trapped 
states become unstable. The horizontal dotted grey line corresponds to E = 0. 



collapse appears first in the self-trapping states. In Fig. 3(a) this corresponds to the 
point where the curve for \Eq\ (dashed line in Fig. 3(a)) terminates at the boundary 
between region I and II (indicated by the vertical dotted line in Fig. 3(a)). At larger 
negative nonlinearities (region I in Fig. 3(a)) the lowest eigenstates invert (as in the 
ID case) and at 5 = the symmetric states also become unstable at #3D = —11.6. In 
Fig. 3(b) we plot typical eigenenergy levels as a function of 5. The curves are similar 
to the ID case except in the limit of large negative nonlinearities (region I), where a 
completely different structure is found. Note that region I cannot be described by the 
two-state model. In this region there is no longer a loop structure as the self trapped 
states are unstable. This parameter region is of interest for investigating tunnelling 
induced collapse, where one begins with a stable symmetric state and adds a potential 
gradient to induce a collapse in one well. However, before discussing the dynamical 
behaviour we consider the stability in the symmetric double-well as a function of the 
barrier height. 
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Figure 4. (a) Critical parameter k cr as a function of h for a spherical geometry 
(A = 1). Stability curves of the symmetric "Jg (curve (ii) in Fig. 1), first excited \& e 
(curve (i)) and asymmetric ground states (iii) are shown as circles, crosses and squares 
respectively, (b) Condensate density 1^(0, z)\ 2 along the z axis as a function of h for 
the state Shown are the central |^(0, 0)| 2 (black) and peak |^(0, z m )\ 2 (grey) 
densities, where z m is the longitudinal position of maximum density in the double-well 
configuration (i.e. centre of each individual well). In both figures, the plotted lines 
connect adjacent data points. The vertical dotted line highlights the critical value h cr , 
above which the stability curve consists of two branches. 



3. Stability of a 3D BEC in a double-well potential 

In this Section we study the stability of a BEC with attractive interactions in a 3D 
symmetric double-well trap as a function of the barrier height h. We solve Eq. (0) 
numerically using the Newton method jHSHHEj- Above the critical atom number A/" cr we 
no longer find stationary solutions. A dimensionless constant k CT relating the scattering 
length a with M CT and the properties of the confining trap, is defined by, [TT] 

— tv cr , yv ) 

a 

where a$ = Jh/muo, m is the mass of the atoms confined in the trap and luq = \ 1 ^ 3 lu± 



is the geometrically averaged trap frequency. 

In Fig. 4(a) we plot the critical constant k cr as a function of h for the case of a 
symmetric double-well trap with A = 1 for the ground and first excited states. We find 
that at a critical value of the barrier height, h cr , there are two branches to the stability 
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curve for the ground state. The upper and lower branches correspond to the symmetric 
(Fig. 1, curve (ii)) and asymmetric (Fig. 1, curve (iii)) eigensolutions respectively, 
(plotted by circles and squares in Fig. 4). Note also that k cr reaches its maximum 
value at a height just above that corresponding to the appearance of the second branch. 
This maximum can be explained by a minimum in the peak density of the double-well 
configuration, which is plotted in Fig. 4(b). As the barrier is raised the condensate 
splits in two, thus reducing its maximum density. However, as the trap splits to form 
two separate condensates, the condensates in each well become compressed and the peak 
density increases again. Our results remain qualitatively unchanged for different values 
of A. This picture is similar to that of Adhikari jHSj, except for the behaviour at large h, 
where Adhikari finds that k CT increases again. As anticipated, we find that, for h > h cr , 
when the system is essentially composed of two separate condensates, the value of k CT in 
each well tends towards the value in a single harmonic trap containing the same number 
of atoms as each half of the double- well. 

Finally, the critical constant k cr as a function of h, for the first excited state, 
\I/ e (crosses in Fig. 4(a)) first increases, reaches a peak value, and then decreases 
approaching the value for the symmetric ground state for large h. This is expected, 
as in the limit of large h, the density distributions of \l/g and \l/ e become very similar. 

4. Tunnelling Dynamics under a time-dependent magnetic field gradient 

The main theme of this Section is to investigate the possibility of observing a tunnelling 
induced collapse. At t < we prepare a stable condensate in a symmetric double-well, 
for a value of the nonlinear constant in region I of Fig. 3(a). In order to induce 
a collapse, a potential gradient is applied at t — 0, i.e., 5 = Rt for t > 0, such that 
the right well has higher potential energy than the left. Subsequently, we study the 
dynamics leading to the onset of collapse, by solving Eq. (Q) numerically. Note that this 
qualitative picture should remain correct, even if 3-body loss terms are included in the 
treatment, although the latter may affect the precise value for the onset of the collapse. 

For attractive interactions, the population difference induced by the addition of the 
gradient does not follow that of the eigenstate, as shown in Fig. 5(a). The effect of 
the nonlinearity is that the ground state is immediately projected onto a superposition 
of states. As the potential gradient is increased, the population in one well reaches a 
critical value and a collapse occurs. The critical gradient corresponding to the collapse, 
shown by the dotted vertical line in Fig. 5, is identified as the point where the interaction 
energy (l-Emtl = 1/2|(7 3D | J |\l/| 4 <ir) diverges, see Fig. 5(b). Note that the critical value 
of the number asymmetry N for which the time-dependent collapse occurs is close to 
the maximum value of |iV max | found for the eigensolution (grey horizontal line in Fig. 
5(a)). 

This prediction becomes clearer if we consider a condensate prepared in the first 
excited state \l/ e , where the time evolution closely follows the eigensolution, see Fig. 
5(b). In this case the collapse occurs at exactly the point where the number asymmetry 
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Figure 5. (a)-(b) Evolution of fractional population difference N as a function of 5 
for a system initially prepared in state (a) ^ g and (b) \I/ e (black lines) . The population 
difference for the eigenstates are also shown as solid grey lines. The vertical dashed 
lines mark the critical gradient at which the eigenstates become unstable (6q = 0.475 
and S' — 0.480 for ty g and ^ e respectively). The vertical dotted lines describe where 
the system collapses in the time-dependent simulation. The collapse occurs when 
\N\ reaches the maximum value -/V ma x| (indicated by horizontal grey lines in (a) and 
(b)) of the number asymmetry predicted by the eigenstates. (c)-(d) Evolution of 
the interaction energy Ei nt (thick black line) when the potential gradient 6 — Rt is 
increased at a constant rate R for (c) ^f g and (d) with corresponding eigenenergies 
shown by grey lines. Other parameters used here: g^r> = —7, h = Ahuj±, A = 1 and 
R = lQ- 3 kujl/a ± . 



exceeds |iV max |, see inset of Fig. 5(b). The critical gradient at which the collapse is 
observed is found to be essentially independent of the rate R at which the gradient is 
increased. 

One can also compare the critical number needed in one well before collapse 
occurs with the prediction for the symmetric potential shown in Fig. 4. By defining 
K CT = A/"4|o|/oo) where M' is the number of atoms in the well which collapses we find 
that K cr = 0.471 and 0.467 for ground and excited states in Fig. 5, which is close to 
the value of k cr = 0.470 predicted by the lower branch of Fig. 4. 

Finally, we discuss typical experimental parameters required for the demonstration 
of the tunnelling induced collapse. In the harmonic oscillator units discussed in 
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Figure 6. Surface plot of the evolution of the density distribution ((^(z, t)\ x 
10 10 cm -3 ) along the z-axis (mm) as a function of time (s) for a BEC initially prepared 
in the ^> g (left) and ^f e (right) states with <?3D = —6. Due to the potential gradient, 
tunnelling is induced to the left (z < 0) well for the ^ g state and to the right well for 
fy e . The condensate instability occurs at t = 0.65 s for \I/ g and 0.67 s for Other 
parameters as in Fig. 5. 



For 7 Li atoms and taking g 3 o = —6 and uj_ = 2tt x 100 Hz, we find M = 1200 which is 
below the critical value for collapse. For an applied field gradient R = (I0~ 3 )(hu±/aj_) 
the collapse occurs at t exp ~ 0.6 s. The collapse is illustrated by the density plots shown 
in Fig. 6. We have confirmed that the collapse can be avoided if the potential gradient 
is ramped up to a value smaller than the critical gradient and then held constant. This 
would not hold if the system were very close to the critical region, in which case number 
fluctuations [SHI ED] could enhance tunnelling and hence induce the collapse at a slightly 
smaller gradient than that predicted by our simple model. 

5. Conclusions 

We have studied the stability of a low temperature atomic BEC with attractive 
interactions in a ID and 3D double-well potential. In particular we highlight a regime 
where the condensate is stable if the population in both wells are approximately equal, 
but becomes unstable if there is sufficient tunnelling from one well to the other. 
We study the dynamics of the system when driven by a time-dependent potential 
gradient and show that a collapse occurs at a critical gradient predicted by the time- 
independent solutions. Although this picture is expected to be qualitatively correct in 
low temperature atomic condensates, further work is required to determine the precise 



Section 2, the number of atoms is given by, 




(7) 
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details of the onset of collapse and the subsequent collapse dynamics, as well as the 
effects of fluctuations [231^01 and finite temperature IHUES- 
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